Exact retrospective Monte Carlo computation of arithmetic average Asian options

Taking advantage of the recent literature on exact simulation algorithms (Beskos et al. [1]) and unbiased estimation of the expectation of certain functional integrals (Wagner [23], Beskos et al. [2] and Fearnhead et al. [6]), we apply an exact simulation based technique for pricing continuous arithmetic average Asian options in the Black & Scholes framework. Unlike existing Monte Carlo methods, we are no longer prone to the discretization bias resulting from the approximation of continuous time processes through discrete sampling. Numerical results of simulation studies are presented and variance reduction problems are considered. Introduction Although the Black & Scholes framework is very simple, it is still a challenging task to efficiently price Asian options. Since we do not know explicitly the distribution of the arithmetic sum of log-normal variables, there is no closed form solution for the price of an Asian option. By the early nineties, many researchers attempted to address this problem and hence different approaches were studied including analytic approximations (see Turnball and Wakeman [20], Vorst [22], Levy [15] and more recently Lord [16]), PDE methods (see Vecer [21], Rogers and Shi [18], Ingersoll [11], Dubois and Lelievre [5]), Laplace transform inversion methods (see Geman and Yor [10], Geman and Eydeland [8]) and, of course, Monte Carlo simulation methods (see Kemna and Vorst [13], Broadie and Glasserman [3], Fu et al. [7]). Monte Carlo simulation can be computationally expensive because of the usual statistical error. Variance reduction techniques are then essential to accelerate the convergence (one of the most efficient techniques is the Kemna&Vorst control variate based on the geometric average). One must also account for the inherent discretization bias resulting from approximating the continuous average of the stock price with a discrete one. It is crucial to choose with care the discretization scheme in order to have an accurate solution (see Lapeyre and Temam [14]). The main contribution of our work is to fully address this last feature by the use, after a suitable change of variables, of an exact simulation method inspired from the recent work of Beskos et al. [1, 2] and Fearnhead et al. [6]. In the first part of the paper, we recall the algorithm introduced by Beskos et al. [1] in order to simulate sample-paths of processes solving one-dimensional stochastic differential equations. By a suitable change of variables, one may suppose that the diffusion coefficient is equal to one. Then, according to the Girsanov theorem, one may deal with the drift coefficient by introducing an exponential martingale weight. Because of the one-dimensional setting, the stochastic integral in this exponential weight is equal to a standard Project team Math Fi, CERMICS, Ecole des Ponts, Paristech, supported by the ANR program ADAP’MC. Postal address : 6-8 av. Blaise Pascal, Cité Descartes, Champs-sur-Marne, 77455 Marne-la-Vallée Cedex 2. E-mails : [email protected] and [email protected]